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The glueball-to-vacuum matrix elements of local gluonic operators in scalar, tensor, and pseu- 
doscalar channels are investigated numerically on several anisotropic lattices with the spatial lattice 
spacing ranging from O.lfm - 0.2fm. These matrix elements are needed to predict the glueball 
branching ratios in J/ip radiative decays which will help identify the glueball states in experiments. 
Two types of improved local gluonic operators are constructed for a self-consistent check and the 
finite volume effects are studied. We find that lattice spacing dependence of our results is very weak 
and the continuum limits are reliably extrapolated, as a result of improvement of the lattice gauge 
action and local operators. We also give updated glueball masses with various quantum numbers. 
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I. INTRODUCTION 

Glucballs, predicted by QCD, are so exotic from the 
point of view of naive quark model that their existence 
will be a direct support of QCD. However, experimen- 
tal efforts in searching for glueballs are confronted with 
the difficulty of identifying glueballs unambiguously, even 
though there are several candidate glueball resonances, 
such as / (1370), /o(1500), / (171Q), and /j(2220), etc.. 
The key problem is that there is little knowledge of the 
nature of glueballs and confined QCD vacuum, which 
requires reliable nonperturbative methods to be imple- 
mented. The numerical study of lattice QCD, which 
starts from the first principles, has been playing an im- 
portant role in this hot field in the last twenty years, 
and extensive numerical studies have been carried out to 
calculate the glueball spectrum Q, H, H, Q . These stud- 
ies give the result that the masses of the lowest-lying 
glueballs range from 1 GeV to 3 GeV, and suggest that 
the J/i/j radiative decays be an ideal hunting ground for 
glueballs. However, apart from the mass spectrum, more 
characteristics are desired in order to determine glue- 
balls in the final states of J/i[> radiative decays, one of 
which is the partial widths of J/if) decaying into glue- 
balls. The first step to estimate these partial widths is 
to calculate the vacuum-to-glueball transition matrix el- 
ements (TME) of local gluonic operators, which are non- 
perturbative quantities and can be investigated by the 
numerical calculation of lattice gauge theory. 

The techniques of lattice calculations in the glue- 
ball sector have been substantially improved in the past 
decade. It is known that very large statistics are nec- 
essary in order for the correlation functions of gluonic 
operators to be measured precisely by Monte-Carlo sim- 



ulation. This prohibits the lattice size being too large. 
On the other hand, because of the large masses of glue- 
balls, the lattice spacing (at least in the temporal di- 
rection) has to be small enough so that reliable signals 
can be measured before they are undermined by statisti- 
cal fluctuations. This dilemma is circumvented by using 
anisotropic lattices, which are spatially coarse and tem- 
porally fine. The potentially large lattice artifacts owing 
to the spatially coarse lattice can be suppressed by the 
implementation of improved lattice gauge actions. These 
techniques has been verified to be efficient in the calcu- 
lations of glueball spectrum and are adopted as the 
basic formalism of this work. 

The glueball matrix elements computed in this work 
are of the form ((% 2 TrG ^ V G pa {x)\G) , where \G) refers to 
the glueball state with specific quantum number J PC — 
++ , 2 ++ , or '", G^(x) is the gauge field strength, and 
g the gauge coupling. The lattice version of the gluonic 
operators <? 2 TrG M „G p(T (x) are constructed by the small- 
est Wilson loops on the lattice. To reduce lattice arti- 
facts, the lattice version of each local gluonic operator is 
improved by eliminating the 0{a 2 s ) [a s here is the spatial 
lattice spacing), while the glueball states \G) are obtained 
by smeared gluonic operators. In order to get a reliable 
continuum extrapolation, five independent calculations 
are carried out on lattices with spatial lattice spacings 
ranging from 0.09 fm to 0.22 fm. As by-products, we 
also calculate the masses of the lowest-lying stationary 
states in all of the symmetry channels allowed on a cubic 
lattice at each lattice spacing. To study the systematic 
error from finite volume, two extra independent studies 
are performed at (3 = 2.4 with the same input parame- 
ters but different lattice size. Note that the local oper- 
ators mentioned above are all bare operators, they need 
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to be renormalized to give the physical matrix elements. 
In this work, the renormalization constants of the scalar 
and tensor gluonic operators at the largest lattice spacing 
are extracted by the calculations of gluonic three-point 
functions on a much larger statistical sample (100,000 
measurements). The pseudoscalar operator renormaliza- 
tion is determined through the calculation of the topolog- 
ical susceptibility. Finally, we give the nonperturbatively 
calculated and phenomenologically significant matrix el- 
ements. 

This paper is organized as follows. In Section II, wc 
give a detailed description of the construction of lattice 
local gluonic operators. Two types of lattice realization 
in each channel are defined, including the improvement 
schemes. The details of the computation, including the 
generation of the gauge configurations, the construction 
of the smeared glueball correlators, the extraction of glue- 
ball masses from the correlation functions, and the deter- 
mination of the lattice spacing, are described in Sec. III. 
In Sec. IV, we estimate the finite volume errors and some 
other systematical errors. The removal of lattice artifacts 
due to the finite lattice spacing, including the extrapola- 
tions to the continuum limit, is described in Sec. V. The 
calculation of three point function and the calculation of 
renormalization constants are described in Sec. VI. Sec. 
VII gives the conclusion and some discussion. 



II. LOCAL GLUONIC OPERATORS AND 
MATRIX ELEMENTS 

In the continuum theory, the lowest-dimensional gauge 
invariant gluonic operators are of the form g 2 TrG M „G pa , 
of which the most commonly studied ones are the scalar 
S(x) (QCD trace anomaly), the pseudoscalar P(x) (topo- 
logical charge density), and the tensor Q^ v (x) (energy- 
momentum density). They all have dimension four and 
are all positive under charge conjugation. The explicit 
forms of them are 

S(x) = g 2 TrG )J , v (x)G^ v (x), 
P(x) = g 2 e IMVprT TrG^(x)G pa (x), 

e„„(a;) - 2g 2r Tr(G^ a {x)G au {x)-U^G 2 {x)). (1) 

If we introduce the chromo-electric and chromo- 
magnetic fields, 

Ei = Gio and Bi = —-tijkGjk, (2) 

the Lorentz scalar S(x) and pseudoscalar P(x) can be 
expressed explicitly by these operators, 

S(x) = 2 3 2 Tr(B(a;) 2 +E(a;) 2 ), 
P(x) = 8 3 2 Tr(E(x) • B(aO), 
Qij(x) = 2g 2 Tr(6fAx) - 6fAx)), (3) 



where 9 and 9 are the traceless, symmetric chromo- 
electric and magnetic tensors, 

9% = TrEiEi - J%TrE 2 

ef s = TrBiBj - i^TrB 2 . (4) 

The scalar operators TrB 2 (x) and TrE 2 (:r), the pseu- 
doscalar TrE(x) • B(x), and the tensors 9 E , 6 B , are all 
irreducible representations of the three-dimensional rota- 
tion group SO(3). 

Denoting the normalized scalar, pseudoscalar, and the 
tensor glueball states as \S), \P), and \Tij), respec- 
tively, the non-zero glueball-to- vacuum matrix elements 
for their annihilation at rest by these operators are 

(27rfS 3 (0)f ( sM) = 9 2 (0\ J d 3 xTrB 2 (x)\S) 

(2^) 3 <5 3 (0)/ (s , B) = g 2 (0\ J d 3 xT r E 2 (x)\S) 

(27r) 3 6 3 (0)f (T . B) = g 2 (0\ J d 3 xTr9 B (x)\T l3 ) 

(27r) 3 <5 3 (0)/ (T , s) = g 2 (0\ J d 3 xT^x)^) 

{2n) 3 5 3 (Q)f (PS) = ,g 2 (0| J d 3 xTrE(x) ■ B(x)\P),(5) 

where no implicit summation is applied. It is straightfor- 
ward to reproduce the matrix elements of the scalar S(x) 
and pseudoscalar operator P{x) by these quantities, but 
for the Lorentz tensor Q pv , the situation is more com- 
plicated. First of all, any hadron state \h) is an eigen- 
state of 0o^, so that the matrix elements (0\®on(%)\h) are 
zero. Together with the traceless property (YJ 6 MM = 0), 

the above condition implies that (0| Qu(x)\h) is zero. 

i 

Therefore, there are only five linearly independent no- 
zero matrix elements of Q^, which can be decomposed 
into the color-magnetic part and the color-electric part, 

mvWTij) = 2(0\9f j (x) 9f j (x)\T ij ). (6) 

Due to the rotational invariance, the five polarizations 
give the same matrix element of the tensor. 

The Lorentz invariance is broken on the spacetime lat- 
tice. The zero-momentum stationary glueball state on 
the lattice must be an irreducible representation (irrcp) 
of the lattice symmetry group, i.e. the 24-element octa- 
hedral point group O. The irreps of group O are classi- 
fied as A\, A 2l E, T\ and T 2 , which are the counterparts 
of angular momentum J for the 5*0(3) rotational group 
and have dimensions 1, 1, 2, 3, and 3, respectively. To- 
gether with parity P and charge conjugate G transforma- 
tions, the full symmetry group of simple cubic lattice is 
an d the total quantum number of the lattice 
glueball state is R PC , where R stands for A\, E, T\ 
or T 2 , and PC can be ++,H — , — h, or . As men- 
tioned above, we are interested in the matrix elements 
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of the scalar, pseudoscalar, and tensor operators, which 
correspond to the irreps R PC through the following rela- 
tion: 0++ <-> 0"+ <-> and 2++ <-> £++®T 2 ++ . 
The zero-momentum operators are obtained by summing 
up all the local operators on the same time-slice in each 
channel. The dimensionless lattice operators are given 
below explicitly. 

The magnetic and electric scalar (labeled (S, B) and 
(S, E) respectively) belong to the Af + representation 
and are defined as 

W*)( t ) = ^£ T rB 2 (x 4 ,i), 



o^ SE \t) 



2„4 



9 a . 
L 3 



^TrE 2 (x,,<), 



(7) 



while the pseudoscalar (labeled (P)) belongs to the A 1 + 
representation and is defined as 

In the above questions, L 3 is the dimensionless spatial 
volume of the lattice and a s the spatial lattice spacing. 
For the spin- two irrep of 50(3), the five polarizations 
are split across the E and T2 irreps of O and so the 
tensor operators must be decomposed into their E and T2 
irreducible contents. The four resulting lattice operators 
are labeled (E,B), (E, E), (T 2 ,B), and (T 2 ,E), and are 
given in terms of their continuum counterparts as 

Of 2 ' B) (*) = 4#E T ^(X^) 



L 3 
L 3 

-,2„4 



(T 2 ,B )(t) = g^Ys^B^t) 



L 3 

M 

L 3 



0^' B \t) = i^^TrEsE^t) 



O, 



(E,E) 



(t) 



L 3 



E 



2V3 



Tt[2£g (x,,t) 



(10) 



where the coefficients guarantee that the five components 
are normalized. Thus we have seven different matrix el- 
ements to be calculated on the lattice, 

T{S,B) = (0|O (S ' B) |A++) 
T(S,E) = (0\O {s ' E) \A++) 



(E,B) ^ (E,B) ]J?++X 



T{E,B) = -(0\O\ 



T(E, E) = ^(0\0{ E ' E) + {E ' E) \E+~ 



T{T 2 ,B) = -<0|Oi 



(T 2 .B) ^ q(T 2 .B) n (T 2 ,B) ]rr++ 



2? 



T{T 2 ,E) = -<0|Oi 
T{PS) = {0\O {PS) \Al 



(T 2 .E) ^ q(T 2 ,E) ^ n (T 2 ,B) |T1+ + 



-+\ 



fll) 



However, in the practical lattice study, the operators 
R (t) ( with i? representing the labels mentioned above) 
are not constructed directly by Bi or Ei, but by the 
proper combinations of different Wilson loops. The first 
definition is based on the linearly combination of a set of 
basic Wilson loops with the requirement that the small- 
a expansion of each combination give the correct contin- 
uum form shown respectively in Eq. (f7jl- l|ll|l . We call the 
local operators through this definition Type-I operators. 
The second construction is that the lattice version of Bi 
and Ei are defined first by Wilson loops and are used 
to compose the local operators according to Eq. (|7l)- (|TTt . 
The local operators through this construction is denoted 
Type-II. 

The following details the constructions of Type-I and 
Type-II local operators on the anisotropic lattice with 
the spatial lattice spacing a s and the temporal lattice 
spacing at- With the tadpole improvement [5|, the tree- 
level Symanzik's improvement scheme is implemented to 
reduce the lattice artifact in defining the local opera- 
tors. Since the aspect ratio of the anisotropic lattice, 
£ = a s /at, is always set to be much larger than 1, the 
leading discretization errors in the local operators are at 
0(a 2 , aj, a s a 2 s ) in this work. 



■,2„4 



0^' E \t) = ^^Tri^x^) (9) 



A. Type-I Operator 



o{ E > B) (t) = ^E^ 2 - B !J(x^) 

Xj 

(E,B ){t) = 9^4 ^ _L Tr j 2B 2 _ B 2 _ S 2 ]( Xi) t) 



The Type-I operators are constructed from a set of ba- 
sic Wilson loops as illustrated in Fig.^ an d Tabled which 
are chosen by the requirement that the first terms of the 
small- a expansion of these loops give the desired contin- 
uum operator forms discussed above. We take following 
steps in the construction. First, this set of Wilson loops 
are acted on by the 24 symmetry operations of the cubic 
point group O (as listed in the Table ITT)) . resulting in 24 
copies with different orientations for each type of Wilson 
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FIG. 1: Wilson loops used in making the Type-I gluonic 
operators. 



TABLE II: The 24 symmetric operations of the cubic point 
group O are represented by the coordinate transformation. 



Index 


Operation 


Index 


Operation 




1 






13 


(x, y, Zj 






y,x) 


2 


(x 7/ z) 


— > ( — Z — 11 — T,\ 


14 


[XV z 




( — a? 


-y, z ) 


3 


(x,y,z) 




15 




-> 


(-y, 


x, z) 


4 


(x,y,z) 


-> {-v,-x,-z) 


16 


(a:,y, 2) 








5 


(x,y,z) 


-> (y,z,x) 


17 


(»,y, z) 




(-X, 


2) y) 


6 


(x,y,z) 


-> (—m,—z,-y) 


18 


(x,y,z) 




{-y, 


—z, x) 


7 


(x,y,z) 


-» (z,y,-x) 


19 


(x,y,z) 




(-X, 


y, -2) 


8 


(x,y,z) 


-* (z,-y,-z) 


20 


(x,y,z) 




(2, - 


y,x) 


9 


(x,y,z) 


-> (y,a;, -z) 


21 


(x,y,z) 




(-2. 


x, —y) 


10 


(x,y,z) 


-> (z, -x, —y) 


22 


(x,y,z) 




(y,- 


x,z) 


11 


(x,y,z) 


-> (a:,z,-y) 


23 


(x,y,z) 




{-y, 


z, —x) 


12 


(x,y,z) 


->• (y> 


24 


(x,y,z) 




{%,- 


z,y) 



TABLE I: The Wilson loops used in the construction of the 
Type-I gluonic operators. They are illustrated by the ordered 
paths in the table, where X,Y,Z, and T are directions and the 
minus sign indicates the path going in the negative direction of 
that axis, (la, It) give the numbers of the spatial and temporal 
links involved in the Wilson loop. These numbers also give 
the powers of the tadpole parameter. 



Index 


Name 


Prototype path 


No. of Links 


i 






(WO 


1 


S-Plaquette 


[X,Y,-X,-Y] 


(4,0) 


2 


S-Rectangle 


[X,X,Y,-X r X,-Y] 


(6,0) 


3 


T-Plaquette 


[X,T,-X,-T] 


(2,2) 


4 


T- Rectangle 


[X J X ) T,-X J -X,-T] 


(4,2) 


5 


S-Chair 


[XY,-X,Z,-Y,-Z] 


(6,0) 


6 


S-Butterfiy 


[X,Y,-X,-Y,Z,-Y,-Z,Y] 


(8,0) 


7 


S-Sunbed 


[X,XY r X-X,Z r Y,-Z] 


(8,0) 


8 


T-Chair 


rx,T,-x,z,-T,.z] 


(4,2) 


9 


T-Sunbed 


[X,X,T,-X,-X,Z,-T,-Z] 


(6,2) 


10 


Knot 


rx j T,-X,-T,Z,-Y,-Z,Y] 


(6,2) 


11 


LS-Knot 


pC,T,-X,-T 1 Z,Z 1 -Y,-Z,-Z,Y] 


(8,2) 


12 


LT-Knot 


[X,X,T,-X,-X,-T,Z,-Y,-Z,Y] 


(8,2) 



loops. Next, they are linearly combined to realize the ir- 
reps of the group O. The coefficients for each irreps, say, 
Ai, E, and T2, are listed in Table UTTI To reduce the lat- 
tice artifact due to the finite lattice spacing, the tree level 
Symanzik's improvement is used, which means that dif- 
ferently shaped Wilson loops are combined to construct 
one operator so that the lattice artifacts are pushed to 
higher order of lattice spacing. Tadpole improvement is 
also implemented to improve the reliability of the lattice 
spacing expansion at the tree level p|. 

Apart from the rotational symmetry, the constructed 
operators should have also definite parity and charge 
conjugation properties. The symmetric/antisymmetric 
combination of a Wilson loop and its parity-transformed 
counterparts gives the positive/negative parity. The 
C = + operators can be realized by taking the real part 
of a Wilson loop. 

Any operator that includes the chromo-electric field 
will involve loops with finite extent in the time direc- 



tion. Operators must be defined for a chosen value of t. 
The simplest way to enforce this definition is to ensure 
the operators on time-slice t are eigenstates of the time 
reversal operators, T about that time slice. The eigen- 
value of this reversal must be the same as the parity of 
the field operator to ensure the correct dimension-four 
operator is reproduced. The chromo-electric scalar and 
tensor operators transform positively under T, while the 
pseudoscalar transforms negatively. The combination co- 
efficients of the time-reversed loops are the products of 
coefficients of original loops and the time reversal eigen- 
values of the time-reversed operator. 

It should be noted that the combination coefficients are 
independent of the shape of loops and each irreducible 
representation corresponds to a specific set of combina- 
tion coefficients, denoted by a[ R \r = 1, 2, . . . , 24. 

For clarity, we give the explicit formula of the thirteen 
zero-momentum gluonic operators as follows: 

(S,B) {t) = ^^■ ReT[[1 _ UciriXtt)] 

rReTr[l - E/ C3 r .0M)]) 
1 

A2u 2 s u t 



r,x 

1 



24u6 

0{ S,E )(t) = £2£ (: ^ ReTr[1 _ C/c3>)t)] 



r,x 

1 



192uf U 2 



ReTr[l - U Ci>r {x,t)\) 



0? E - B) (t) = (MX^ReTVf/c^M) 
T ReTiU C2tr {x,t)) 



r,x 

1 

96u 6 ' 



< (2,*)^ReTW c ., |P (a:,t)} 



r,x 

1 



96u^ ReTlUc ^ t)) 
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TABLE III: Combinational coefficients used in Eq. 112H to construct the irreps of the cubic point group. 



reO 




af + (l,l) 


a? + (l,2) 


«f + (2,l) 


af + (2,2) 


a r 2 (2,1) 


rj-,+ 

a r 2 (2,2) 


a r 2 (2,3) 


1 


1 


-1 


i 

73 


i 

75 


1 





1 





2 


1 


-1 


1 

73 


i 

1 

75 


-1 





1 





3 


1 


1 


1 

75 


-1 








1 


4 


1 


1 


1 

f 3 


1 

73 


1 








1 


5 


1 





7/3 


2 

7"7 





1 








6 


1 





2 

75 


2 

75 





1 








7 


1 


-1 


1 

75 


1 

75 


-1 





-1 





8 


1 


-1 


1 

75 


1 

75 


1 





-1 





9 


1 


1 


i 

75 


i 

7S 


1 








-1 


10 


1 


1 


1 

75 
2 

73 




-1 








-1 


11 


1 





2 

73 





-1 








12 


1 





2 
~i 


2 
77! 





-1 








13 


1 


-1 


1 

75 


I 

73 


-1 





-1 





14 


1 


-1 


1 

75 


1 

73 


1 





-1 





15 


1 


1 


1 

75 


1 

73 
1 

73 


1 








-1 


16 


1 


1 


i 

75 


-1 








-1 


17 


1 





2 

75 


2 

73 





-1 








18 


1 





2 

75 


2 

73 





-1 










1 


1 

-1 




1 

J5 


1 

1 


n 

u 


-1 

1 


u 


20 


1 


-1 






-1 





1 





21 


1 


1 






-1 








1 








73 


73 






22 


1 


1 






1 








1 


23 


1 












1 








24 


1 





75 


73 





1 









o 



(Tt 2 ,b: 



(*) 



0f T2 ' E \t) 



0< ps >(t) 



E^ 2 M(- 

r,x 



11 

96^ 



96u[ 

e 2 E 



:ReTrUc e r (x,t) 



ar 2 (2,t)( 



96u*u 2 



■ReTrt/cr fx t) ^ rotation on the i-th prototype loop in Table [J The 

tadpole parameters u s and u t here come from the renor- 
malization of spatial and temporal gauge links, respec- 
ReTrC/ C7 ,.(a;,i))tively: Uj(x) -> Uj(x)/u s and U t (x) -> U t (x)/u t . The 
determination of it s and «t is described in Section III. 
B. Type-II Operators 



1 



48uf 
ReTr(7 G 



Xx,t) 



;^ ! ^ReTrt/ Cl0 , r (x,t) 



96ufttf 



1 



ReTrU Cllir .(x,t) 
jReTrU Cl2 Jx,t)), 



192u s s u 2 

where Uc ir (x,t) is the loop generated by operating 



(12) 



Generally speaking, there can be many ways to define 
the lattice local operators as long as these definitions give 
the same continuum limit, and some may have smaller 
lattice artifacts at finite lattice spacing compared to oth- 
ers. This motivates us to design another type of lattice 
gluonic operators, called Type-II operators in this work, 
apart from the construction described above. Both types 
of operators are used for a self-consistent check. 

According to the non-Abelian Stokes theorem [l^, a 
rectangle Wilson loop P"* b (x) of size a x b, with a,b 
small, can be expanded as 



■f + <tl> ( F^(x) + ~(a£>„ + bD v )Fp,{x) + ^(2a 2 ^ + 2,abD„D v + 2b 2 D 2 v )F„ v {x) 



— (a 3 Dl + 2a 2 bDlD v + 2ab 2 D^D 2 v + b 3 D 3 v )F^(x] 



6 



{abf 



+ -F pv {x){aD p + bD v )F pv {x) + —F^(x)(a 2 Dl + b 2 D 2 u )F^(x) 



0((ab) 4 ). 



(13) 



For simplicity, the factor ig is absorbed into the quantity 
F pv and will be reconsidered when comparing with the 
continuum form. This expression can be simplified by 
the clover-type combination which is defined as 



C<£; b) (x) 



1 



(P$ b (x) + P^(x) + Plf_ v {x) + Plll{x)) 
+ (P^ a (x) + P b v ™{x) + P b _™ v {x) + P h _X{x)) 
] (14) 



where 



P±n±v{x) = U± ll {x)U± v {x ± afx)U^Jx ± bv)U\. v (x). 

(15) 

The small ab expansion of P±n±v(x) is similar to Eq. I|13(l 
by replacing a and b with ±et and ±6, respectively. This 
clover-type combination is illustrated in Fig. O 



J 



(a) C 



(i.D 



(b) C 



,(1,2) 



FIG. 2: The clover-shape combinations of spatial plaque- 
tte and rectangle, which are used to derive the gauge field 
strength F M „ on the lattice. Here i and j are the indices of 
the spatial direction. 

As a result, the tree- level expansion of C pu (x) is ex- 



plicitly derived as 



1 



ReC M „(x) = l+-(abYFUx) 



1 



+ -^{abYia' + b z )F^{x){D'+ Di)F^ 



+ 0((ab)% 
ImC^x) = abF^ v {x) 



+ -ab(a 2 + b i )(D^ + Di)F^{x) 
+ h a b) 3 (F flv ) 3 + .... 



(16) 

Using lmCui/(x) with different a and b as the ele- 
mentary components, the local operators F^ v F pa can 
be defined through the lattice version of the gauge field 
strength tensor 



F^ix) = -(5ImCtf(x)-ImCtf(x)) 

= a 2 (f mv (x) - ^a\F^(x)) 3 + 0(a 6 )) , 

(17) 

which is improved up to 0{a A ) by the combination of 
the clover-leaf diagram in Fig. EI a) and the wind-mill 
diagram in Fig- EJb) . Thus the fourth-dimensional gauge 
operator F^ v F p(7 can be derived from F pl/ (x)F pcr (x), 

F liV {x)F pa {x) = a i F pv (x)F pa (x) 



: a 8 F pu (F 2 „(x) + F 2 pa (x)) F pa (a 



+ 



(18) 



In order to check the effectiveness of the improvement 
scheme of lattice local operators, we have also derived 
another lattice definition of F pv F pa (x), which is through 
the direct combination of Im C^ v (x) , say, 



F, v (x)F pa (x) = 1 (71m Ci' 1 (x) ■ Im C# (x) - Im (x) ■ Im C 1 * (x) - Im C 1 * (x) ■ Im C% (x)) 

= a^F^F^x) - \a*F p „ (F 2 v (x) + F 2 a {x)) F pa (x) - ^{D 2 p + D 2 v )F pv {D 2 p + D 2 a )F pa + 



(19) 



I 

One can find that, at tree level, the lowest order difference between these two definitions, denoted by A, is 

A = -^ a8 ( D l + D 2 u )F pu (D 2 p + Dl)F pa + .... (20) 
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In the practical calculation, based on these two ap- 
proaches, we construct two version of operators in 
irreps. In Section IV, we find the different definitions re- 
sult in about 3-4% discrepancy for the measured matrix 
elements. 

The discussion above are based on the classical (or tree 
level) series expansion of Wilson loops. For this to be re- 
liable, the tadpole improvement should be applied, which 
means the tadpole parameter should be included in the 
above expressions. Specifically, a rt-link spatial Wilson 
loop C^ v (x) should be divided by a tadpole factor 
From the data analysis to be shown in Sec. |H] that the 
tadpole improvement alleviates most of the dependence 
on finite lattice spacing. The improvement scheme of the 
local operators corresponding to EiEj is a little differ- 
ent from that of BiBj. Since at <C a s , all the tempo- 
ral Wilson loops included in the improvement have only 
one lattice spacing extension in the time direction. In 
other words, for the temporal loops, we do not include 
the windmill diagrams which involve two lattice spacing 
in the time direction. Thus the combination coefficients 
in the above expression are modified accordingly. We 
omit the explicit expression here. 



III. Numerical Details 

Since the implementation of the tadpole improved 
gauge action on anisotropic lattices was verified to be 
very successful and efficient in the determination of the 
glueball spectrum Q, we use the same techniques to 
calculate the glueball matrix elements. We adopt the 
anisotropic gauge action used by Morningstar and Pear- 
don in 



s 



5 £l sp 4 £f2 tp 



i a 



sir 



3 uivft 



12 £u% 12 u\u\ 



s 



}, 

(21) 

where (3 = Q/g 2 , g is the QCD coupling constant, £ is 
the aspect ratio for anisotropy, u s and ut are the tadpole 
improvement parameters of spatial and temporal gauge 
links, respectively, and Oc's are the sums of various Wil- 
son loops over the total lattice (the explicit expression 
can be found in Ref. In practice, u s is defined by 

the expectation value of the spatial plaquette, namely, 



u s = ((l/STriV)) 1 / 4 , and u t is set to 1. Theoreti- 
cally, the bare anisotropy £ should be finely tuned to 
give the correct physical anisotropy £p%s = a s /at, but in 
our practical case £ is always taken as the same as t; p hys 
because the discrepancy is shown to be within 1 — 2 per- 
cents when the improved action in Eq. H21(l is used Q. 
For each coupling constant f3 and £, u s is determined 
self-consistently in the Monte Carlo updating. 

The gauge configurations were generated by using 
Cabibbo-Marinari(CM) pseudo-heatbath and the SU(2) 
subgroup micro-canonical over-relaxation (OR) methods. 
Three compound sweeps were performed between mea- 
surements, where a compound sweep is made up of one 

TABLE IV: The input parameters for the calculation. Values 
for the coupling (3, anisotropy £, the tadpole parameter iQ, the 
single-link smearing parameter X s , the double-link smearing 
parameter A/, lattice size, and the number of measurements 
are listed. 
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i 


ui 




A/ 


L 3 


x T 




2.4 


5 


0.409 


0.1 


0.5 


8 3 


x 40 


999 x 100 












12 3 


x 64 


100 x 100 












16 3 


x 80 


100 x 100 


2.6 


5 


0.438 


0.1 


0.5 


12 3 


x 64 


86 x 100 


2.7 


5 


0.451 


0.1 


0.5 


12 3 


x 64 


100 x 100 


3.0 


3 


0.500 


0.4 


0.5 


16 3 


x 48 


100 x 100 


3.2 


3 


0.521 


0.4 


0.5 


24 3 


x 72 


79 x 100 



CM updating sweep followed by 5 OR sweeps. The mea- 
surements of n m b configurations are averaged in each bin, 
and ribin bins are obtained. Table IIVI lists the relevant 
input parameters for lattices with 5 different (3. For the 
case of (3 — 2.4, there are three lattice volumes to study 
the finite volume effects. In order to calculate the ma- 
trix element such as (0\O\G), it is desirable to have the 
glueball state \G) determined as precisely as possible. In 
this work, the glueball states \G) with quantum number 
R = Af c , A% c , E PC , Tf c , and T 2 PC are generated by 
smeared gluonic operators Og , which are constructed by 
exploiting link-smearing and variational techniques in a 
sequence of three steps outlined below. First, for each 
generated gauge configuration, we perform six smear- 
ing/fuzzing schemes to the spatial links, which are vari- 
ous combinations of the single-link procedure 



Uj(x) - P su(3) {U 3 (x) + X s U k (x)Uj(x + k)Ul(x + j)}, 

±(Mj) 

Uj(x) = Psu<${U3{x)Uj(x + 3) + \f U k (x)U j (x + k)U j (x + j + k)U k (x + 2j)}, (22) 

±(fe?y) 

i 



where Psu{3) denotes the projection into SU (3) and is realized by the Jacobi method [8| in this work. 
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7 S 



r 




FIG. 3: Wilson loops used in making the smeared glueball 
operators. 



The six schemes are given explicitly as s\ 



h t <8> s a 8 > A/ <8> s \ s > A/ ® s a 3 > where s/f denotes 
smearing/fuzzing procedure defined in Eq. (|22[) . A s and 
Xf are tunable parameters for smearing and fuzzing and 
take the optimal value A s = 0.1 ( 0.4 at £ = 3) and 
Xf = 0.5. Secondly, for each smearing/fuzzing scheme, 
we use ten Wilson loops illustrated in Fig. which are 
the same as used in |4| , as prototypes to construct the op- 
erator <f>^ (t) which is a linear combination of different 
oriented spatial loops, invariant under spatial transfor- 
mation, and transforms according to the irreps R. There 
are four independent constructions for each R (except 
) for each smearing/fuzzing scheme, thus the 



for 



operator Og(t) is a linear combination of 24 operators, 



O 



(R)AR) 



(t). The coefficients 



Va are deter- 



( R) 

mined by a variational method so that O s projects 
mostly to a specific glueball states \G). 

What we obtain from the MC simulation are the 24 x 24 
correlation matrix 

C^(t)=^(0|^)(t + r)4 fl) (r)|0), (23) 



where the vacuum subtraction 

^(*) = 4 fl) (*)-<o|4 fl) (*)|o) 



(24) 



is only applied to the channel which has a vacuum 
expectation value. The coefficients are determined in the 
data analysis stage by minimizing the effective mass 

rh{t D ) = -—In— — m m z , (25) 



a/3 



where the time separation for optimization is fixed to 
tr> = 1. This is equivalent to solving the generalized 
eigenvalue equation 



TABLE V: The lattice spacings determined by r^ 1 = 410(20) 
MeV at different p. 



p 


ro/a s 


a s /ro 


a s (fm) 


2.4 


2.17(1) 


0.461(2) 


0.222(1) 


2.6 


2.74(2) 


0.365(2) 


0.176(1) 


2.7 


3.09(2) 


0.326(2) 


0.156(1) 


3.0 


4.05(4) 


0.247(2) 


0.119(1) 


3.2 


4.76(5) 


0.210(2) 


0.101(1) 



and the eigenvector Vq R \ corresponding to the lowest 
effective mass fh^to), yields the coefficients v K 0a ' for the 
operator (t) which, under ordinary circumstances, 
best overlaps with the lowest-lying glueball Go in the 
channel of interest. The operators most overlapping with 
the excited states can also be constructed accordingly. 

After the smeared operators are determined, the glue- 
ball masses and the matrix elements we are concerned 
with can be extracted by fitting two two-point functions, 
namely, the smeared-smeared correlation function, 



C'ss(t) 



J2(0\Os(t + r)Os(r)\0) 

T 

E^^<°l°^>HOs|o) e --' ' 



1 



2m G V 

and the smeared-local one 



|(G|O s |0)| 2 e- MG * (t-oo), (27) 



C SL {t) = ^(0|O s (i + r)O L (r)|0) 

T 



-m n t 



l 



2M G V 



(G|O s |0)(0|O i |G)e- MG * (i-oo). 



(28) 



Css(t) and CsL(t) can be fitted simultaneously with 
three parameters, i.e., (0\Os\G), (0\Ol\G), and the 
ground state glueball mass Mq. (0|Ol|G) is the glue- 
ball matrix element we would like to obtain. Of course, 
before the physical results can be derived, a proper renor- 
malization scheme should be performed. 

In the following, we describe the calculation details 
step by step. 



A. Setting the Scale Using the Static Potential 

The five lattice spacings a s are determined by calcu- 
lating the heavy quark static potential. This part of the 
calculation is independent of the production runs. The 
static-quark potential V(r) can be extracted from the 
averages of Wilson loops W(r, t), 



C{t D )w 



(R) 



-t D m(t D ) 



G(0)v 



(i?) 



(26) 



W(t, t) = Z(r) exp(-*V(r)) + 



(29) 
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which can be measured precisely on the lattice. For 
each /3, 200 configurations are generated, each of which 
is separated by 30 compound sweeps so that the auto- 
correlation effects are reduced. Secondly, these config- 
uration are fixed to temporal gauge so that the Wilson 
loops can be calculated more easily. Different smearing 
schemes are applied to Wilson loops with different size to 
increase the overlapping with the ground states. Finally, 
the static potential is fitted by the model 



TABLE VI: The fitted ground state masses in the twenty 



V(r) = V& + OT + — , 

r 



(30) 



with three parameters Vo, <J, and e c in the correlated fit 
method. The lattice spacing a s is determined by 



ro 



1.65 + e c 



(31) 



according to the relation r 2 ' dV '(r)/ 'dr] r=ro = 1.65, where 
ro is hadronic scale parameter jg. The lattice spacings 
for different (3 are listed in Table E| 



B. The Glueball Mass Spectrum 

Glueball masses can be obtained by fitting the two- 
point functions Css{t) directly. It is taken as a testi- 
mony of the effectiveness of the improvement and smear- 
ing scheme 4J. Even though we are interested in the 
glueball matrix elements of the four channels , , 
E ++ , and in this work, glueball masses of all the 20 
channels R PC are calculated as a by-product. 

After the implementation of variational-optimization 
in each channel, we can obtain a specific operator which 
overlaps most with the a specific state and has little con- 
taminations from other states with the same J PC , so that 
we can use a single mass term 



n !+\ 7 ( -Mt i -M(T-t)\ 

Css{t) = Z ( e + e >) 



(32) 



to 



fit the two-point function in a time range 
i , . . . , t max , which can be determined by observing 



the effective mass plateau. As a convention in this work, 
we use Mg to represent the mass of a glueball state in 
the physical units and M to represent the dimensionlcss 
mass parameter in the data processing with the relation 
M = Meat- Generally speaking, for most channels in 
each of the five (3 cases, the overlap of the specific oper- 
ators on the ground states are all larger than 90%. 

Before performing the continuum extrapolation, we 
check the finite volume effects of glueball masses. Three 
independent calculations at (3 — 2.4, £ = 5 were carried 
out on a 8 3 x 40 lattice, a 12 3 x 64, and a 16 3 x 80 lat- 
tice. These lattices have spatial volumes of (1.76 fm) 3 , 
(2.64 fm) 3 , and (3.52 fm) 3 , respectively, with the lattice 
spacing a s ~ 0.22 fm from r^ 1 = 410(20) MeV. For these 
three runs, all the input parameters are the same except 
the different lattice volume. 



R channels at j3 
L 3 x T = 16 



2.4 and £ = 5 on three different lattices 



3 x 80, 12 3 x 64, and 



x 40. 





L = 16 


L = 12 


L = 8 




0.308(3) 


0.308(2) 


0.312(2) 




1.082(14) 


1.051(13) 


1.079(10) 


. |_ 

A + 


0.604(11) 


0.618(4) 


0.618(3) 




1.03(3) 


1.075(16) 


1.01(3) 




0.825(8) 


0.818(8) 


0.820(5) 




0.805(6) 


0.804(8) 


0.808(5) 


A 1- 


1.047(10) 


1.067(12) 


1.053(12) 




0.993(11) 


0.976(10) 


0.994(27) 


E++ 


0.542(3) 


0.536(2) 


0.541(2) 


E+- 


0.919(17) 


0.957(8) 


0.960(6) 


n h 

E + 


0.699(4) 


0.698(4) 


0.695(3) 


E~~ 


0.879(5) 


0.891(7) 


0.882(5) 


T++ 


0.826(5) 


0.832(4) 


0.834(3) 


T+- 


0.657(6) 


0.661(5) 


0.663(4) 


Tf+ 


0.932(5) 


0.940(5) 


0.936(4) 


r— 


0.865(12) 


0.884(5) 


0.893(4) 


rp + + 

J 2 


0.542(3) 


0.536(3) 


0.538(2) 


T +- 


0.807(4) 


0.808(11) 


0.799(8) 


T - + 


0.697(5) 


0.700(4) 


0.700(2) 


T 2 — 


0.903(4) 


0.893(6) 


0.896(3) 



The calculated glueball mass spectra are listed in Ta- 
ble El where one can find at a glance that the finite vol- 
ume effects are very small and in most case the changes 
are within errors and consistent with zero statistically. 
More precisely, we use the following scheme to illustrate 
the finite volume effect (FVE) quantitatively. Let M de- 
note the average value of the glueball masses from the 
three lattice volumes, M(L) denotes the glueball mass 
measured on lattice L 3 x T. The fractional change of 
the glueball mass is defined by S G (L) = 1 - M(L)/M. 
The results for these fractional changes are shown in 
Fig. 0] Each point in the figure shows the fractional 
change 5g{L) of the glueball mass, the error bars come 
from the statistical errors of M{L), the solid lines indi- 
cate 8 a = 0, the dotted lines above and below the solid 
lines indicate Sq = 0.02 and 5q = —0.02, respectively. 
All changes are statistically consistent with zero, sug- 
gesting that systematic errors in these results from finite 
volume arc no larger than the statistical errors for these 
physical volumes. Since the physical volumes of the other 
(3 values are not smaller than the 8 3 x 40 at (3 — 2.4, we 
shall neglect the FVE for the higher (3 results. 

The fitted glueball masses at different coupling con- 
stant (3 are listed in Table IVIII where the statistical er- 
rors are also quoted. The dimensionless products of r 
and the glueball masses Mq are shown as functions of 
(a s /ro) 2 . To remove discretization errors from our esti- 
mates, the results for each level in these figures must be 
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1.. 



I x X nz_ _ 



16 12 8 16 12 8 16 12 8 
++ +- -+ 



16 12 8 



FIG. 4: Finite-volume effects of glueball masses at /3 = 2.4, 
£ = 5. Each point shows the fractional change Sg(L) = 1 — 
M(L)/M of the glueball mass, where M(L) is the glueball 
mass measured on the lattice L 3 x T with L — 8, 12, and 16, 
and M is the average values over those from different lattices. 
The errorbars come from the statistical errors of M(L). The 
lattice labels L are shown along the horizontal axis, and the 
labels along the vertical axis are the irreps of lattice symmetry 
group. The solid lines indicates 5g = 0, and the dotted lines 
above the solid lines indicates So = 0.02, and the dotted line 
lines below the solid lines indicate Sg = —0.02. 



extrapolated to the continuum limit a s /ro — > 0. From 
perturbation theory, the leading discretization errors are 
expected to be 0(a 2 , a^, a s a 2 s ). As discussed in Ref. Q, 
the 0(a 2 ,a s a 2 ) errors could be negligible compared to 
the 0(a A s ) errors for most calculated glueball masses ex- 
cept for Af + glueball, so the fit model for these glueball 
masses is chosen to be 



M G (a s )r = Af G (0)r o + c 4 ^f 



(33) 



The masses versus (a s / r o) 2 i as wei l as the fitted curves 
are shown in Fig. El and [S] One can find from 
the figures that the data obey this function very well. 
However, for Af + glueball mass, the 0(a 2 s ) error seems 
still very large. So we keep the linear term of a 2 in the 
fit model, namely, 



M G (a s )r = M G (0)r + c 2 ^ + c 4 - 



(34) 
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0.05 0.1 0.15 0.2 0.25 0.3 

(a s /r ) 2 

FIG. 5: Masses of PC = H — h glucballs in terms of ro against 
the lattice spacing square (a s /ro) 2 . The calculated values are 
plotted in points, and the curves are the best fits. 
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FIG. 6: Masses of PC = — h glueballs in terms of ro against 
the lattice spacing square (a s /ro) 2 . The calculated values are 
plotted in points, and the curves are the best fits. 



with (3 = 2.4, 2.6, 2.7, 3.0 and 3.2. For comparison, we 
have listed the results from the previous work |4|. 



Already there has been some discussion on the possible 
reason for this large discretization error in scalar channel; 
one can refer to Ref. Q for details. 

We list several continuum extrapolated glueball masses 
in Table IVIIII We note that the earlier work Q was 
carried out with = 1.7, 1.9, 2.2, 2, 4, 2.5 and 3.0. In the 
present work, we concentrate on finer lattice spacings 



C. The Glueball Matrix Elements 

As described above, we use the correlated minimal- 
X 2 method to fit the smeared-smeared and smeared-local 
two-point functions, namely, Css(t) an d Cshit), simul- 
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FIG. 7: Masses of PC = H — glueballs in terms of ro against 
the lattice spacing square (a s /ro) 2 . The calculated values are 
plotted in points, and the curves are the best fits. 
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FIG. 8: Masses of PC = glueballs in terms of ro against 

the lattice spacing square (a s /ro) 2 . The calculated values are 
plotted in points, and the curves are the best fits. 



taneously using the fit functions 
Css(t) 

C SL (t) = XY (e- Mt + e- M( - T ~t) 



= X 2 e -Mt + e -M(T-t) 



(35) 



where X = {0|O g (0)|G) / \J2M G V and Y = 
(0\O L (0)\G)/V2M G V. Here V is the spatial lat- 
tice volume in physical units, say, V — L 3 a^, and 
Mq = Mat . As an example, we show the effective 
mass plot at (3 = 2.4 (L = 8 lattice) in Figs. |Hl and 



TABLE VII: The fitted ground state masses M in the twenty 
R PC channels at five different /3'b. The errors quoted here are 
all statistical errors from correlated minimal-x 2 method. 



R 



PC 



-at 



At + 



E+^ 
E+- 

E— 

r+- 



rp + + 

J 2 

T +- 

T -+ 
1 2 

T — 



= 2.4 /3 = 2.6 (3 = 2.1 j3 = 3.0 (3 = 3.2 



0.312(2) 
1.079(10) 
0.618(3) 
1.01(3) 

0.820(5) 
0.808(5) 
1.053(12) 
0.994(27) 

0.541(2) 
0.960(6) 
0.695(3) 
0.882(5) 

0.834(3) 
0.663(4) 
0.936(4) 
0.893(4) 

0.538(2) 
0.799(8) 
0.700(2) 
0.896(3) 



0.264(2) 
0.851(17) 
0.470(4) 
0.897(20) 

0.665(4) 
0.644(7) 
0.850(7) 
0.748(5) 

0.418(2) 
0.745(7) 
0.543(4) 
0.695(6) 

0.645(5) 
0.526(6) 
0.737(7) 
0.685(5) 

0.419(2) 
0.631(4) 
0.540(3) 
0.705(5) 



0.247(2) 
0.741(11) 
0.410(3) 
0.815(5) 

0.593(6) 
0.560(5) 
0.745(26) 
0.672(9) 

0.374(3) 
0.663(13) 
0.486(3) 
0.623(10) 

0.561(3) 
0.459(6) 
0.664(9) 
0.610(7) 

0.375(4) 
0.568(7) 
0.479(4) 
0.642(5) 



0.325(3) 
0.961(10) 
0.515(6) 
1.036(12) 

0.779(17) 
0.721(11) 
1.02(4) 
0.827(21) 

0.475(5) 
0.843(13) 
0.627(3) 
0.797(29) 

0.736(8) 
0.596(9) 
0.844(11) 
0.776(8) 

0.476(8) 
0.720(7) 
0.604(10) 
0.805(23) 



0.279(5) 
0.822(13) 
0.439(4) 
0.879(19) 

0.683(17) 
0.616(8) 
0.868(18) 
0.721(9) 

0.414(4) 
0.734(16) 
0.514(3) 
0.682(13) 

0.628(4) 
0.510(5) 
0.714(6) 
0.650(9) 

0.413(5) 
0.617(4) 
0.510(3) 
0.680(10) 



1101 where the smeared-local correlation functions are 
plotted on the negative-time side of the t axis. In each 
channel, the effective mass plateau for Css are flatter 
than that of Csl- Nevertheless, both of them can be 
fitted by the same mass parameter consistently. For all 
the five /3's, the fit parameters X,Y, effective masses 
M, as well as the fit windows, are listed in Table IIXI 
H EU Km and IXlTTl In these tables, the first group 
of data is for Type-I operators, the second group for 
Type-II operators, and the last two row of data for the 
second definition (Eq. I|19l) ) of Type-II operators in T2 
irreps. We give a brief interpretation of the meaning 
of the value of X parameter here. During the practical 
performance of the variational method discussed in Sec. 
Ill, the eigenvectors v( R > in Eq. l|26|l are normalized as 

v' R ) C(0)v( R ' — 1, so that the parameter X indicates 
the relative overlap of the smeared operator to the 
specific stationary state (the ground state here). In 
these tables, all the X parameters are very close to 
1 which implies that the smeared operators couple 
almost completely to the ground state and there is little 
contamination from excited states. The x 2 's per degree 
of freedom of the data fitting are all acceptable and also 
listed in the tables. 

From the measured matrix elements of the two defini- 
tions of the Type-II local operators in tensor channel (T 2 
irreps), we can get an estimate of the 0(a A s ) lattice arti- 
facts. Recalling the discussion in Sec. II, after restoring 
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TABLE VIII: Continuum-limit glueball masses Ma. The cor- 
responding results in Reference 4| are also quoted for com- 
parison. 



R PC 


Possible J PC 


voMg (this work) 


roMa [4] 


Af + 


0++ 


4.16(11) 


4.21(11) 


E++ 


2++ 


5.82(5) 


5.85(2) 


r r+-\- 


2++ 


5.83(4) 


5.85(2) 




3 ++ 


9.00(8) 


8.99(4) 


T++ 


3 ++ 


8.87(8) 


8.99(4) 




0"+ 


6.25(6) 


6.33(7) 


T+~ 


1+- 


7.27(4) 


7.18(3) 


E-+ 


2" + 


7.49(7) 


7.55(3) 


rn—+ 


2~ + 


7.34(11) 


7.55(3) 


T+- 


3+" 


8.80(3) 


8.66(4) 




3+" 


8.78(5) 


8.66(3) 


T~ 


\— 


9.34(4) 


9.50(4) 


E~~ 


2~~ 


9.71(3) 


9.59(4) 


T 2 - 


2~ 


9.83(8) 


9.59(4) 




3~" 


10.25(4) 


10.06(21) 


E+- 


2+" 


10.32(7) 


10.10(7) 


At~ 


o+- 


11.66(7) 


11.57(12) 
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FIG. 9: The effective mass plateaus of the smeared-smeared 
and smeared-local correlation functions at (3 — 2.4, £ = 5 on 
the lattice 8 3 x 40. The smeared-local correlation functions 
are plotted on the negative-time side of the t axis. The Type-I 
local operators are used. 

the fact ig to F^, the difference of the two definition of 
the operators (T2,B) is 

A = ^ 8 s9 2 (Dl + Dl)F^(D 2 p + Dl)F pa + ..., (36) 

with the indices /x, v, p, and a varying accordingly. The 
ratio of the matrix elements of the two definitions of op- 
erator (T2,B) is plotted in Fig. with respect to the 
lattice spacing. It is known that this difference comes 
totally from the lattice discretization and will disappear 
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FIG. 10: The effective mass plateaus of the smeared-smeared 
and smeared-local correlation functions at /3 = 2.4, ( = 5 on 
the lattice 8 3 x 40. The smeared-local correlation functions are 
plotted on the negative-time side of the t axis. The Type-II 
local operators are used. 



in the continuum limit. This is the case from the figure. 
Even though the relative difference is less than 3-4%, 
the deviation of the ratio from 1 is detectable on coarse 
lattice with larger lattice spacing and decreases when ap- 
proaching to the continuum limit. On the finest lattice 
we are using, the difference is consistent with zero within 
the error. This result shows that, as far as the Type- 
II operators are concerned, after the implementation of 
Symanzik's improvement scheme along with the tadpole 
improvement, the measured matrix elements of the two 
versions of lattice local operator have a few percent dif- 
ferences at finite a s , but approaches the same continuum 
limit as a s goes to zero. 

Here comes the discussion of the FVE of the matrix 
elements. By analogy with the FVE analysis of glueball 
masses, the glueball matrix elements are also calculated 
at = 2.4 and £ = 5 on a 8 3 x 40 lattice, a 12 3 x 64, and 
a 16 3 x 80 lattice. For these three runs, all the input pa- 
rameters are the same except the different lattice volume. 
The extracted matrix elements are listed in Table IXlVl 
IX VI and IX VII respectively for Type-I and Type-II opera- 
tors. As illustrated in the tables, the finite volume effects 
of matrix elements are very small and in most case the 
changes are consistent with zero statistically. We also de- 
fine the relative deviation Sa(L) = 1 — /g(L)//g to show 
the FVE quantitatively, where fo is the mean value of the 
matrix elements of glueball G averaged over the three lat- 
tice volumes, fa{L) denote the matrix element of glueball 
G measured on lattice L 3 x T. The results of Sc{L) are 
shown in Fig. ^| In the figure, each point shows Sq{L) 
of the matrix element of a glueball G, and the error bars 
come from the statistical errors of fo(L). The solid lines 
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TABLE IX: The matrix elements extracted at j3 = 2.4 (on lattice 8 s x 40). The smeared-smeared(SS) and smeared-local(SL) 
correlators are fitted simultaneously in the time window t(SS) and t(SL). The fitted masses M, X, and Y parameters are listed. 
The first group of data is for Type-I operators, the second group for Type-II operators, and the last two row of data for the 
second definition (Eq. 1191 ) of Type-II operators in T2 irreps. 



R 


t(SS) 


t(SL) 


M 


Y( x 10" 2 ) 


X 


X 2 /d.o.f 


(S,B) 


1-4 


1 - 3 


0.311(2) 


21.3(1) 


0.994(3) 


0.92 


(S,E) 


1 - 4 


1-4 


0.311(2) 


21.1(2) 


0.994(3) 


0.69 


(PS) 


1 - 5 


2-4 


0.617(4) 


5.24(11) 


0.985(4) 


1.59 


(E,B) 


1 - 7 


1 - 5 


0.541(2) 


2.64(3) 


0.995(2) 


0.64 


(E,E) 


1 - 7 


1-4 


0.541(2) 


1.52(9) 


0.995(2) 


0.64 


(T 2 ,B) 


2-7 


1-4 


0.539(3) 


1.84(12) 


0.989(2) 


0.57 


(T 2 ,E) 


2-7 


1 - 3 


0.539(3) 


0.85(5) 


0.990(2) 


0.52 


(S,B) 


1 - 9 


1 - 7 


0.311(2) 


11.54(6) 


0.993(3) 


1.12 


(S,E) 


1 - 9 


1 - 7 


0.311(2) 


16.38(13) 


0.995(3) 


0.75 


(PS) 


2-7 


1 - 6 


0.617(4) 


5.33(11) 


0.985(4) 


1.36 


(E,B) 


1 - 7 


1 - 6 


0.541(2) 


2.50(2) 


0.995(2) 


0.89 


(E,E) 


1 - 7 


1 - 3 


0.542(2) 


1.35(4) 


0.995(2) 


0.87 


(T2,B) 


2-7 


1 - 5 


0.541(2) 


2.24(1) 


0.990(3) 


1.34 


(T2,E) 


2-7 


2-4 


0.542(2) 


0.73(4) 


0.989(2) 


0.45 


(T2,B) 


2-7 


2-5 


0.543(2) 


2.37(1) 


0.990(2) 


1.06 


(T2,E) 


2-7 


1-4 


0.542(2) 


0.73(4) 


0.989(2) 


0.46 



TABLE X: The same as in Table but for /3 = 2.6 



R 


t(SS) 


t(SL) 


M 


Y( x 10 -2 ) 


X 


X 2 /d-o.f 


(S,B) 


1 - 9 


1 - 7 


0.264(2) 


6.48(5) 


0.994(3) 


0.92 


(S,E) 


1 - 9 


1 - 7 


0.264(2) 


6.37(11) 


0.994(3) 


0.69 


(PS) 


2-7 


1 - 5 


0.470(4) 


2.14(3) 


0.985(4) 


1.59 


(E,B) 


1 - 7 


1 - 6 


0.424(2) 


0.716(10) 


0.995(2) 


0.64 


(E,E) 


1 - 7 


1 - 3 


0.423(2) 


0.39(4) 


0.995(2) 


0.64 


(T2,B) 


2-7 


1 - 5 


0.421(1) 


0.562(7) 


0.989(2) 


0.57 


(T 2 ,E) 


2-7 


1-4 


0.421(2) 


0.23(2) 


0.990(2) 


0.52 


(S,B) 


1 - 9 


1 - 7 


0.265(2) 


4.01(3) 


0.993(3) 


1.12 


(S,E) 


1 - 9 


1 - 7 


0.263(2) 


5.06(3) 


0.995(3) 


0.75 


(PS) 


2-7 


1 - 6 


0.470(5) 


2.16(3) 


0.985(4) 


1.36 


(E,B) 


1 - 7 


1 - 6 


0.424(2) 


0.706(13) 


0.995(2) 


0.89 


(E,E) 


1 - 7 


1 - 3 


0.424(2) 


0.40(3) 


0.995(2) 


0.87 


(T2,B) 


2-7 


1 - 5 


0.423(2) 


0.698(5) 


0.990(3) 


1.34 


(T 2 ,E) 


2-7 


2-4 


0.421(2) 


0.22(2) 


0.989(2) 


0.45 


(T2,B) 


2-7 


2-5 


0.422(2) 


0.719(10) 


0.990(2) 


1.06 


(T2,E) 


2-7 


1-4 


0.421(2) 


0.22(2) 


0.989(2) 


0.46 



indicates So = 0, the dotted lines above and below the 
solid lines indicate 5q — ±0.02. The largest effects from 
the finite volume appear to occur in the TYPE I opera- 
tor of (E,B). All deviations are statistically consistent 
with zero, suggesting that systematic errors in these re- 
sults from finite volume are no larger than the statistical 
errors. Since the physical volumes of the other (3 values 
are not smaller than the 8 3 x 40 at (3 — 2.4, we shall 
neglect the finite volume effects for the higher (3 results. 

Based on the discussions above, in the following contin- 
uum extrapolation of matrix elements, we use the results 
calculated on the lattice 8 3 x 40 at (3 = 2.4. As for the 
Type-II operators, we take the results with the lattice 
operators defined by F^ l/ (x)F pa (x) 

Using the fitted M and Y parameters listed in Ta- 
ble [K|2^n] as well as the physical values of the lattice 



spacings a s listed in Table IVI for each (3, the lattice ma- 
trix elements (0\Ol\G) can be obtained by (0\Ol\G) — 
\[TMqVY . As a result, the matrix elements in units of 
Tq 3 are listed in Table IXVlII and illustrated by FigurelTSI 
It is clear that the continuum limits can be extrapo- 
lated neither by the function with a single a 2 s term nor 
with a single af term. We fit them with the form 

T(a s )rl = T(0)rl + c 2 ^ + c^, (37) 

where T is the final non-renormalized continuum limit 
results of the glueball matrix elements. 

Note that these matrix elements are all calculated by 
bare lattice operators. Before the renormalization of the 
local operators is performed, we cannot draw any con- 
clusions of physical interest at present stage. However, 



14 



TABLE XI: The same as in Table [DO but for /3 = 2.7 



R 


i(SS) 


t(SL) 


M 


Y( x 10" 2 ) 


X 


X 2 /d.o.f 


(S,B) 


2-9 


1 - 7 


0.247(2) 


5.03(4) 


0.990(3) 


0.92 


(S,E) 


2-9 


1 - 7 


0.246(2) 


4.88(13) 


0.991(3) 


0.50 


(PS) 


1 - 7 


2-6 


0.431(3) 


1.88(5) 


0.991(2) 


1.07 


(E,B) 


3-8 


1-4 


0.376(3) 


0.575(10) 


0.982(4) 


0.75 


(E, E) 


3-8 


1 - 3 


0.376(4) 


0.27(3) 


0.981(4) 


0.39 


(T 2 ,B) 


4-8 


2-6 


0.379(5) 


0.423(10) 


0.977(7) 


1.14 


(T2,E) 


4-8 


1-4 


0.378(6) 


0.20(2) 


0.975(8) 


1.39 


(S,B) 


2-9 


3-7 


0.246(2) 


3.21(4) 


0.990(3) 


0.86 


(S,E) 


2-9 


3-7 


0.247(2) 


3.80(10) 


0.990(3) 


0.58 


(PS) 


2-7 


1 - 6 


0.431(2) 


1.90(5) 


0.991(2) 


1.08 


(E,B) 


3-8 


2-5 


0.376(3) 


0.592(10) 


0.981(4) 


0.60 


(E,E) 


3-8 


2-5 


0.376(3) 


0.27(2) 


0.982(4) 


0.22 


(T2,B) 


4-8 


2-6 


0.386(5) 


0.546(10) 


0.988(8) 


1.79 


(T 2 ,E) 


4-8 


1 - 3 


0.378(5) 


0.19(2) 


0.975(8) 


1.28 


(T2,B) 


4-8 


3-6 


0.386(5) 


0.565(9) 


0.987(8) 


1.82 


(T 2 ,E) 


4-8 


2-3 


0.378(5) 


0.19(2) 


0.975(8) 


1.29 






TADTf VTT. 

IKdLiSLi All: 


The same as 


m lable I1AI but lor p — o.O 






R 


t(SS) 


t(SL) 


M 


Y( x 10" 2 ) 


X 


X 2 /d.o.f 


(S,B) 


2-9 


1-4 


0.324(3) 


1.56(3) 


0.985(3) 


0.99 


(S,E) 


2-9 


1-4 


0.326(4) 


1.50(4) 


0.985(3) 


0.82 


(PS) 


2-5 


1-4 


0.542(6) 


0.706(15) 


0.981(6) 


2.00 


(E, B) 


3-7 


1-4 


0.480(6) 


0.176(7) 


0.967(8) 


0.22 


(E,E) 


3-7 


1 - 3 


0.480(6) 


0.107(14) 


0.967(8) 


0.25 


(T 2 ,B) 


3-8 


1 - 3 


0.487(5) 


0.137(5) 


0.969(7) 


1.28 


(T 2 ,E) 


3-8 


1 - 3 


0.485(6) 


0.060(9) 


0.967(7) 


0.94 


(S,B) 


2-9 


2-7 


0.326(4) 


1.07(2) 


0.984(3) 


0.68 


(S,E) 


2-9 


1 - 8 


0.324(3) 


1.09(3) 


0.986(3) 


0.49 


(PS) 


2-5 


1-4 


0.541(6) 


0.715(15) 


0.981(5) 


1.84 


(E,B) 


3-7 


1-4 


0.477(6) 


0.179(5) 


0.964(8) 


0.72 


(E,E) 


3-7 


1 - 3 


0.480(6) 


0.086(5) 


0.967(8) 


0.22 


(T 2 ,B) 


3-8 


1-4 


0.490(5) 


0.172(4) 


0.971(7) 


1.55 


(T2,E) 


3-8 


1 - 3 


0.485(6) 


0.053(6) 


0.968(8) 


1.17 


(T 2 ,B) 


3-8 


1-4 


0.489(5) 


0.179(4) 


0.971(7) 


1.46 


(T 2 ,E) 


3-8 


1-3 


0.485(6) 


0.053(6) 


0.968(8) 


1.17 








1 he same as 


m laoiG iiA-i out tor p — o.z 






R 


t(SS) 


t(SL) 


M 


Y( x 10 -2 ) 


X 


X 2 /d.o.f 


(S,B) 


4-7 


1 - 6 


0.279(6) 


0.518(11) 


0.966(9) 


1.74 


(S,E) 


4-7 


1-4 


0.278(5) 


0.534(26) 


0.965(9) 


0.18 


(PS) 


3-9 


3-6 


0.439(8) 


0.294(15) 


0.932(9) 


0.38 


(E, B) 


3-7 


1-4 


0.417(4) 


0.062(4) 


0.967(3) 


2.77 


(E, E) 


3-7 


1 - 3 


0.415(4) 


0.048(8) 


0.965(3) 


1.49 


(T 2 ,B) 


3-9 


1-4 


0.414(3) 


0.0548(13) 


0.959(4) 


0.67 


(T 2 ,E) 


3-9 


1 - 3 


0.415(3) 


0.033(5) 


0.960(5) 


0.67 


(S,B) 


4-7 


1 - 8 


0.282(5) 


0.391(7) 


0.968(7) 


1.98 


(S,E) 


4-7 


1 - 6 


0.280(5) 


0.376(11) 


0.967(8) 


0.31 


(PS) 


3-9 


3-6 


0.439(8) 


0.299(15) 


0.932(9) 


0.40 


(E,B) 


3-7 


1 - 3 


0.415(4) 


0.0693(23) 


0.966(4) 


3.29 


(E,E) 


3-7 


1 - 3 


0.416(4) 


0.037(4) 


0.965(4) 


3.05 


(T 2 ,B) 


3-9 


1 - 7 


0.416(4) 


0.068(2) 


0.961(6) 


0.96 


(T 2l E) 


3-9 


1 - 3 


0.415(3) 


0.028(3) 


0.960(4) 


0.56 


{T 2 ,B) 


3 - 9 


2 - 5 


0.415(4) 


0.068(3) 


0.960(5) 


0.74 


{T 2 ,B) 


3 - 9 


1 - 3 


0.415(3) 


0.028(3) 


0.960(4) 


0.54 
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TABLE XIV: The finite volume effects of matrix elements at 
(3 = 2.4 and £ = 5 on three different lattices L 3 xT — 16 3 x 80, 
12 3 x 64, and 8 3 x 40. The fitted Y parameters are listed in 
normal numbers. For clarity, the Y parameters for lattice 
L = 12 and L = 8 are re-scaled to L — 16 due to the relation 
Y = (0\Ol{0)\G}/^/2M g V and illustrated by bold numbers. 
This table lists the Y parameters of Type-I operators. 





L — 16 


L = 12 


L = 8 


(S,B) 


7.53(5) 


7.50(5) 


7.52(5) 








91 97f 1 fil 


(S,E) 


7.38(12) 


7.51(11) 


7.47(8) 




7 91 


11 C >fi/'171 


91 19(931 


(E,B) 


0.956(12) 


0.993(8) 


0.936(11) 






1 ^(191 


9 R4(3l 


(E,E) 


0.56(5) 


0.52(4) 


0.54(3) 




0.56(5) 


0.81(6) 


1.52(9) 


(T2,B) 


0.633(8) 


0.649(8) 


0.651(4) 




0.633(8) 


1.000(12) 


1.841(11) 


(T2,E) 


0.29(2) 


0.27(3) 


0.30(2) 




0.29(2) 


0.42(4) 


0.85(6) 


(PS) 


1.85(4) 


1.88(4) 


1.85(4) 




1.85(4) 


2.89(6) 


5.24(11) 



DC 




0.05 0.1 0.15 0.2 0.25 0.3 



(a s /r o) 



FIG. 11: The data points represent the ratios of the matrix 
elements of the two definitions of operator (T2,B) at different 
lattice spacings. The straight line is plotted to show how far 
the ratios deviate from 1. 



TABLE XVI: The same as Table IXTVl but for the second 
definition of Type-II operators Q(T 2 ,B) and 0{T 2 ,E). 





L = 16 


L = 12 


L = 8 


(T 2 ,B) 


0.825(6) 


0.835(8) 


0.837(4) 




0.825(6) 


1.285(12) 


2.367(12) 


(T 2 ,E) 


0.26(2) 


0.26(2) 


0.26(2) 




0.26(2) 


0.40(3) 


0.73(4) 



TABLE XV: The same as Table IxWI but for Type-II oper- 
ators. 





L = 16 


L = 12 


L = 8 


(S,B) 


0.410(3) 


0.408(3) 


0.408(2) 




4.10(3) 


6.29(4) 


11.54(6) 


(S,E) 


5.73(6) 


5.82(4) 


0.579(4) 




5.73(6) 


8.97(6) 


16.4(1) 


(E,B) 


0.885(11) 


0.895(9) 


0.883(7) 




0.885(11) 


1.378(14) 


2.50(2) 


(E,E) 


47(2) 


0.47(2) 


0.48(2) 




0.47(2) 


0.72(3) 


1.35(4) 


(T2,B) 


0.783(5) 


0.792(4) 


0.790(4) 




0.783(5) 


1.220(6) 


2.24(1) 


(T 2 ,E) 


0.25(2) 


0.26(2) 


0.26(2) 




0.25(2) 


0.40(3) 


0.73(4) 


(PS) 


1.87(4) 


1.90(4) 


1.89(4) 




1.87(4) 


2.92(6) 


5.33(11) 



we can give some comments on the different behaviors 
of the matrix elements of the Type-I and Type-II oper- 
ators from the continuum extrapolation. Based on our 
experience from the calculation of glueball masses, the 
continuum SO(3) symmetry is approximately restored 
for all the lattice spacings we use in this work, since the 
calculated glueball masses in + and E ++ channel are 
coincident. The left panel of Fig.^|is the plot of the ma- 
trix elements of the Type-I operators, where one can find 
that the calculated glueball matrix elements of T2 and 
E irreps do not show this symmetry restoration. This is 
probably due to the definition of the Type-I operator in- 
troduced in Sec. II. The Type-I operators with different 
quantum number are defined by the real part of different 
Wilson loops composed of different numbers of spatial 
gauge links, such that the overall tadpole improvement 
factors (different powers of the tadpole parameter u s ) are 
different for T2 and E representation. The conjectured 
power counting of the tadpole parameter u s for Wilson 
loops is a naive approximation and may bring additional 
deviation to the local operators. We have carried out 
the test that the naive tadpole parameters are replaced 
by the vacuum expectation value (VEV) of the corre- 
sponding Wilson loops. As we expected, the two matrix 
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TYPE-I TYPE-II(1) TYPE-II(2) 

(ps) ;;;; n;;; it ;;;; ;;r;; ;;;;; ;;;; ;;;; i ;;;; gg ;;;; ;;;; ;;;; ;;;;; ;;;; ;;;; ;;;; ;;;;; 



(T,E) 



FIG. 12: Finite-volume effects on the results of the (3 — 2.4, 
£ = 5 calculation. Each point shows the fractional change 
Sg(L) = 1 — /g(L)//g in the matrix element of a glueball 
G, where /g(L) is the matrix element of G measured on the 
lattice L 3 x T with L = 8, 12, and 16, and fa is the average 
values over those from different lattices. The errorbars come 
from the statistical errors of fc(L). The lattice labels L are 
shown along the horizontal axis, and the labels of different 
local operators are specified along the vertical axis. The solid 
lines indicates Sg = 0, and the dotted lines above the solid 
lines indicates 8g = 0.02, and the dotted line lines below 
the solid lines indicate S G = -0.02. The FVE of different 
definitions of local operators are also shown in the figure. 



elements tend to coincide better. In contrast with the 
Type-I operators, all the Type-II operators are defined 
by the improved lattice version of the gauge strength 
tensor, F^, and thus have the same correction factor 
coming from the tadpole improvement. The right panel 
of Fig. 1131 shows the behaviors of the matrix elements of 
the Type-II operators. It is clearly seen that the approxi- 
mate SO(3) symmetry restoration takes place for T 2 and 
E representation. 

The matrix elements of phenomcnological interest, de- 
noted s, p and t, are defined as 



(2tt) 3 <5(0)s 
(27r) 3 5(0)p 



(27r) 3 5(0)te^ = 



d 3 xS(x)\G), 
d 3 xP{x)\G), 
d 3 x& flu (x)\G), (38) 



where e M „ = e ufi , e w = 0, e^e^ = 1, and S(x), P{x), 
(x) are defined in Eq. and Eq. J3J). Recalling 
that (27r) 3 <5(0) is replaced by L 3 a 3 on finite lattices, and 
combining Eq. I|7I11|I . s, p, and t can be reproduced by 
the calculated matrix elements listed in Table IX VIII as 



TABLE XVII: Final results of matrix elements of lattice glu- 
onic operators. The data in table are in the unit of r^ 3 with 
r^ 1 = 0.410(20) MeV. 



(3 = 2.4 



(3 = 2.6 
(L = 12) 



13 = 2.7 
(L = 12) 



13 = 3.0 
(L = 16) 



13 = 3.2 
(L = 24) 



T{S, B)i 
1 (b, E)i 


87(2) 
oo(z ) 


90(2) 

CO/ A \ 

oy(4j 


94(3) 

9t(0 ) 


92(3) 
yu(o ) 


85(4) 

on/ ff\ 

88(7) 


T(E,B)i 
1 (E, E)i 


14.2(3) 
8.2(0) 


12.5(5) 
6.8(tU) 


13.2(7) 


12.6(8) 

*7 f\{ 1 0\ 
i ) 


12.8(10) 
9.6(22) 


T{T 2 ,B)! 

1 {l2,E)i 


9.9(3) 

A K f A \ 

4.0(4J 


9.8(4) 

A 1 { P.\ 

4.t(0) 


9.7(5) 

A fcfQ\ 

4.o(oJ 


9.8(7) 

A O ( T\ 

4.3(7) 


11.0(6) 
O. /(14) 


T{PS)i 


30(1) 


40(2) 


46(3) 


52(3) 


61(5) 


T{S,B) n 
T{S,E)n 


47(1) 
67(2) 


56(1) 
70(2) 


60(2) 
71(4) 


63(2) 
65(4) 


64(3) 
62(4) 


T(E,B)u 
T{E,E)u 


13.4(4) 
7.3(3) 


12.3(4) 
7.0(7) 


13.6(5) 
6.2(6) 


12.8(8) 
6.2(6) 


13.8(9) 
7.4(10) 


T{T 2 ,B) n 
T(T 2 ,E)u 


12.7(3) 
3.9(3) 


12.6(5) 
3.8(5) 


13.0(7) 
4.4(5) 


12.8(7) 
3.8(6) 


13.6(9) 
5.6(10) 


T(PS)h 


31(1) 


40(2) 


47(2) 


53(3) 


62(5) 



Wo 
*r 2 rg 



8T(PS), 

2\T(E,B)-T(E,E)\, 
2\T(T 2 ,B)-T(T 2 ,E)\. 



(39) 



Fig. [H and Table IXVIIII illustrate the behaviors of 
these matrix elements with respect to the lattice spac- 
ings. We only show the data reproduced from the Type- 
II operators, which can be only renormalized in this work 
(see Section V). The continuum limits of tx 2 and tE are 
consistent within error bars. The deviation of the central 
values comes mainly from the discrepancy of T(E, E) and 
T(T 2 ,E). 



V. The Non-Perturbative Renormalization of Local 
Gluonic Operators 

The lattice local gluonic operator O and its continuum 
counterpart O con t is related by the renormalization con- 
stant Zo{a), 



Z o (a)0. 



(40) 



2(T(S,B)+T(S, E)), 



The key question is to choose a proper normalization con- 
dition, so that the renormalization constant can be de- 
termined. In this work, we choose the energy-momentum 
tensor in the glueball state as the normalization condi- 
tion. 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.05 0.1 0.15 0.2 0.25 0.3 0.35 

(a s /r ) 2 (a s /r ) 2 

FIG. 13: Glueball matrix elements(ME) in terms of rj) against the lattice spacing square (a s /ro) 2 . The calculated values of 
ME's are plotted in points, with different colors corresponding to different channels. The curves are best fits to the calculated 
results by the model T(a B )rQ = T(0)rj) + C2(a s /ro) 2 + Cn(a s /ro) i . The left panel illustrates the result of Type-I operators, while 
the right panel is for Type-II operators. 



TABLE XVIII: The non-renormalized matrix elements s, p, and t in units of r^ 3 are reproduced from the lattice results of 
type-II operators at different lattice spacings. The errors quoted here are purely statistical. The continuum limits are also 
listed for each channel. 





P = 2.4 


P = 2.6 


P = 2.7 


P = 3.0 


P = 3.2 


Continuum 


sr'o 


228(5) 


252(5) 


262(9) 


256(9) 


252(10) 


227(7) 


pro 


248(8) 


320(16) 


376(16) 


424(24) 


496(40) 


589(43) 


t E 4 


12.2(5) 


10.6(8) 


14.8(8) 


13.2(1.0) 


12.8(1.4) 


13.6(4.1) 


tx 2 r z 


17.6(0.8) 


17.6(0.8) 


17.2(0.9) 


18.0(1.0) 


16.0(1.4) 


15.8(1.9) 



The energy-momentum tensor for the pure gauge the- 
ory 

= = --S^TrF 2 + 2TrF„ a F va (41) 

satisfies d^T^ — and does not need an overall renor- 
malization in the continuum. More specifically, Too = 
Tr(E 2 - B 2 )). At the classical level, is trace- 

less. When quantum corrections are included, the renor- 
malized energy-momentum tensor, T Ml/ , takes a non- 
vanishing trace part, T^ v , which comes from the anoma- 
lous breaking of scale invariance and is called the QCD 
trace anomaly, 

^ = ~4^T^ F2 = ^I^rr(. 9 2 B 2 + 5 2 E 2 ), 

9 ^ (42) 

where (3(g) = - P g 3 / (4tt) 2 is the (3 function of QCD to 
the lowest order of the coupling constant g with /3q = 
11 in the pure gauge case. Thus, T^ v can be written 
explicitly as the sum of the traceless and trace parts, 

T M „ = % v + f„ v . (43) 
On the other hand, T^ v defines the Hamiltonian oper- 



ator of the theory 

H = J ' d 3 xT 00 (x,0), (44) 

which is finite and scale independent. In the rest frame 
of a glueball, the matrix element of the Hamiltonian H 
in the glueball state is just the glueball's rest mass, 

(G(M G , p = 0)| / d 3 zT o(x, 0)\G(M G , p = 0)) 
G (G\G) 

(45) 

If the glueball states \G(E,p)) are normalized as 

(G(E p ,p)\G(E p ,p')) = 2E p (2ir) 3 5 3 (p-p'), (46) 

we have (G(M G ,p = 0)| J d 3 xT 00 (x, 0)\G(M G ,p = 0)) = 
2M g (2it) 3 5 3 (0) = 2p°p°{2ir) 3 5 3 (0), which implies that 

(G(p)\ J d 3 xT^0)\G( P ))=2p^(2ir) 3 S 3 (O). (47) 

Combining Eq. I|45l) and l|47[) . we get the normalization 
conditions in the glueball's rest frame, 
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(G(p)\ I d 3 xf^(x,0)\G(p)) = 
(G(p)\ J d 3 xf^(x,0)\G(p)) = 

According to the definitions in Eq. El and Eq. the 

E and T 2 operators are related to Ty by 

{E,E) _ (E,B) k f n -T 22 , 
(E,E) _ (E.B) k 2?33 _ j,^ _ 
(T 2 , S )_ (T 2 ,B) ^ £ .. fcf . fcj (4Q) 

whose matrix elements in a zero-momentum glueball 
state vanish, thus we cannot renormalize E ++ and T% 
operators directly by the normalization condition in 
Eq.0S]in the glueball rest frame. A possible way around 
this difficulty is to assume the almost restoration of 
Lorentz invariance at the lattice spacings we use in this 
work, so the overall renormalization constant of the lat- 
tice version of can be determined by one of its com- 
ponents. In fact, this is justified to some extent by two 
facts. First, it is argued that the rotational invariance can 
be restored if the scale parameter z = m A ++L is larger 



2(pfp- - -M G S^)(27r) 3 S 3 (0), 
^M G 5^(2n) 3 S 3 (0). 

(48) 

I 

than 10 [T^j, an d our smallest lattice gives the value 
z = 0.31 x 5 x 8 = 12.4 which meets this requirement. 
Secondly, in our lattice calculations of the mass spectrum 
and matrix elements, the coincidence of E ++ channel and 
channel implies that this rotational restoration is ac- 
tually realized. Based on the discussion above, we choose 
the component Too to do the renormalization of the ten- 
sor operator in this work. In the practical study, we 
calculate the matrix elements (G|0±|G)iat, where 0± is 
the operator g 2 Tr(E 2 ±B 2 ) and \G) = \A+ + ), \E ++ ), or 
|T 2 f+ ). As we have addressed before, 0_ is proportional 
to Too, and + is proportional to the trace anomaly To . 
Thus, the renormalization constants of the scalar and 
tensor operators can be extracted from these matrix ele- 
ments. 

The matrix elements (G\0±\G)i a t can be obtained by 
calculating the three-point function 



C 3 (ti,t2) - (0lOa(-ti)O±(0)O.(fe)|0) = X) zp f T/2 

m,n 4£/ « £/m V 

1 



(0\O s (-h)\n){n\O±(0)\m)(m\O s (t2)\0) 



= £ 



4E n E m V 2 



\O s (0)\n)(n\O ± (0)\m)(m\O s (0)\0)e 



-E n ti e -E m t 2 



1 



{2M G vy 



\O s (0)\G)\ 2 (G\O ± (0)\G)e- MG(tl+t ^ (h, t 2 -> 00) 



r 



(50) 



where Os (t) is the smeared zero- momentum operator simultaneously with the fitting model 
which generates the ground state |G) from the vacuum, 
and Mq the ground state mass. Using the asymptotic 
form of the two-point function 



v 2 -Mt 

X e 



Css(t) 

C 3 (t,t ) = X 2 We- M(t+t0 \ 



C'ss(t) 



1 



(53) 



2MqV 



\O s (0)\G)\ 2 e 



2-M G t 



(t — > 00) (51) where to = 1. IT is related to the matrix element 



(G\O ± (0)\G) by 



and dividing the three-point function by proper two-point 
function, one gets 



(G\O±(0)\G) = 2M G VW. 



(54) 



1 



C 3 (t u t 2 ) _ 

C S s{h) ~ 2M G V 



(G\O ± (0)\G)e- MGt2 (ti,t 2 ->oo). 

(52) 

In the practical data processing, by analogy with the ex- 
traction of the matrix element, we fit C2(t) and C 3 {t, to) 



From Eq. ijiDjl. Eq. and the fitted matrix elements 
(G\0±\G)\ a t, the renormalization constants for scalar 
and tensor gluonic operators can be derived as 



1 



11 



Zs(a s ) 



32ir 2 M 2 aj 



(G\0 + \G) 



lat 1 
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TABLE XIX: The matrix elements of operator 0± (Type-II) 
in glueball states are listed. S represents the scalar glueball 
state, while T(E) and T(T2) denotes the tensor glueballs in 
E and T2 irreps, respectively. The renormalization constants 
extracted from these matrix elements are given also. 



p(a s )/p(0) 



\G) 


(G\0+\G)u t 


Zs 


<G|0-|G)l at 


Zt 


S 


32(9) 


1.1(3) 


13(5) 


0.7(3) 


T{E) 


102(16) 


1.0(2) 


51(15) 


0.53(15) 


T(T 2 ) 


101(16) 


1.0(2) 


53(15) 


0.51(15) 



Z T (a s ) 



9 2 *M%ai 



(G\0-\G) 



lai 1 



(55) 



where the coupling constant g 2 comes from the relation 
f 00 = (l/g 2 )O-. 

Unfortunately, the gluonic three-point function is far 
more noisy than the gluonic two-point function in Monte 
Carlo calculation. In this work, the renormalization of 
gluonic operators is performed only at j3 = 2.4 on the 
lattice 8 3 x 40 with £ = 5. As many as 100,000 measure- 
ments are carried out, the signals of three point functions 
are still weak with large fluctuation. In the practical 
computation, the matrix elements of Type-I and Type- 
II operators of 0( S ' S ) and 0^ S,E ' are all calculated. It 
is found that the three-point functions involving Type-I 
operators are so noisy that the matrix elements cannot 
be extracted reliably. The three-point functions involv- 
ing Type-II operators behave better, from which we ob- 
tain the renormalization constants. The matrix elements 
(G|0±|G)iat and the resultant renormalization constants 
at = 2.4 are listed in Table IXIXl 




0.25 



(a s /r ) 



FIG. 14: The non-renormalized matrix elements s, p and t 
from type-II operators are plotted versus lattice spacing. The 
continuum extrapolations are also carried out. 

As for the renormalization of the pseudoscalar opera- 
tor, we tentatively use the phenomenological value of the 




0.25 



(a s /r r 



FIG. 15: The a s -dependences of Xi an d p(fls) are plot- 



ted. For comparison, x^ 2 (a 3 ) and p(a s ) are rescaled by their 
continuum-extrapolated values, say, ~x^ 2 (Q) andp(0), respec- 
tively. One can find that the a s -dependences of Xz/ 2 ( ffls ) and 
pro(a s ) are very similar and their ratios are fairly constant. 



topological susceptibility, x 1 / 4 = 180 MeV, as the nor- 
malization condition. The quantity x is defined by 



X = d 4 x{q(x)q(0)), 



(56) 



where the topological charge density q(x) is proportional 
to the pseudoscalar operator P(x) (defined in Eq. (JJJ 
as q(x) = -^jP(x). The lattice version of x, denoted 
by Xl i i s defined by the lattice pseudoscalar operator 
0(PS)(x) as 



XL 



(32tt 2 ) 2 L*Ta 4 s 



-£(0( ps \x)0( ps \y)), (57) 



where L and T are the lattice sizes in spatial and tem- 
poral directions, respectively. 

We have calculated xl(°s) for all the five (3, at each 
f3 20000 measurements are carried out. In Table IXXl are 
listed the results of X// 4 (°s) in units of MeV for differ- 
ent /3, as well as the non-renormalized continuum value 
Xz/ 4 (0) after the continuum extrapolation a s — * 0. It is 
obvious that the value of Xi/ 4 (°s) increases along with 
the decreasing of lattice spacing. For comparison, the 
non-renormalized matrix elements pr^ (a s ) at different /3's 
are also listed in Table IXXl In Fig. ^| are plotted the 
a s -dependences of x}/ 2 (fls) and precis) (rescaled by their 
continuum-extrapolated values, respectively), as well as 
their ratios, where one can find that the a s -dependences 
(a s ) and pr^ia^ are very similar and their ratios 
are fairly constant. Using x 1 ^ 4 = 180 MeV, the renor- 
malization constant Zp of pseudoscalar operator in the 
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TABLE XX: The non-renormalized topological susceptibility Xl{o-s) calculated at different lattice spacing are listed in physical 
units. Also listed are the ratios of x^ 2 {a s ) and the non-renormalized matrix elements of pseudoscalar p(a s ) after rescaled by 
their continuum extrapolation values Xi/ 2 (0) and p(0). 



p 


2.4 


2.6 


2.7 


3.0 


3.2 


continuum 


xi /4 (a s )(MeV) 

p(a s ) x ' /2 (0) 
X^ /2 (a s )p(0) 


242(4) 
248(8) 


277(5) 
320(16) 


299(5) 
376(16) 


323(6) 
424(24) 


351(8) 
496(40) 


391(15) 
589(43) 


1.10(7) 


1.08(9) 


1.09(8) 


1.06(10) 


1.05(11) 


1.00(13) 



continuum limit can be extracted as 



Z P = 



,1/2 



xi /2 (o) 



0.21(2) 



(58) 



VI. RESULTS AND DISCUSSION 

After the continuum extrapolation and the renormal- 
ization of the local operators, we can discuss the physical 
implication of our lattice results. 

As we addressed above, we have tried to extract the 
renormalization constants for scalar and tensor operator 
only at (3 = 2.4. In fact, the signals of the three-point 
functions of type-I operator are very poor, so the renor- 
malization constants Z$ and Zt are obtained only for 
Type-II operators. We are lucky with this situation be- 
cause all the Type-II operators are made up of the lattice 
version of the gauge field strength, F^, thus they all have 
the same normalization constant, say, the same tadpole- 
improvement factor. Therefore, for Type-II operator, the 
Zt extracted from T 00 can be applied to other compo- 
nents involved in the glueball matrix elements calculated 
in this work. 

Considering Eq. (JSJ and © , the non-renormalized ma- 
trix elements of the S(x), P{x), and 0(x) in the con- 
tinuum can be obtained by the lattice results, and are 
shown in Table IX Villi We notice that, although the 
(S, B) and (S, E) matrix elements have sizeable a s de- 
pendence for the Type-II operator, as seen in Fig. 1131 
the total scalar matrix element, which is the sum of the 
two, is much flatter in (a s /ro) 2 . This is also the case for 
the tensor matrix element. With the observation that 
the non-renormalized matrix elements of the scalar and 
the tensor depend on the lattice spacing (a s /ro) 2 very 
mildly, we speculate that there is not a large lattice spac- 
ing dependence in the renormalization constant, and will 
use Z$ and Zt computed at (3 = 2.4 as an approxima- 
tion of the renormalization constants in the continuum 
limit. We shall check this in the future when computer 
resources are available for high statistics calculation at 
larger (3. 

Taking the average of the renormalization constant Z$ 



from Table IXIXI we get a continuum extrapolated value 
for the matrix element (0|5|0++) = 15.6 ± 3.2 GeV 3 . 
Based on the scaling properties of QCD and trace 
anomaly, both the QCD sum rule and the soft meson 



theorem lead to an estimate that relates the scalar 
glueball matrix element to the gluon condensate, 



<0|S|0++) =16n\^M G , 



(59) 



where Go = (Q| ^f-G£„G" \0) is the gluon condensate, 
b = (ll/3)iV c -(2/3)JV/ and M G the scalar glueball mass. 
Taking N f = 0, G = 0.012 GeV 4 , and M G w 1.7 GeV, 
this matrix element is estimated to be ~ 6.3 GeV, which 
is about two and a half times smaller than our lattice 
result. This discrepancy might be attributable to the 
fact that the quenched lattice calculation gives a gluon 
condensate which is about 0.14 ± 0.02GeV 4 [l^. This is 
larger by an order of magnitude than that used in QCD 
sum rule. If the relation Eq. I|59|) still holds in the pure 
gauge theory, using the quenched gluon condensate, the 
estimated scalar matrix element is estimate to be (21 ± 
l)GeV 3 which is in good agreement with our quenched 
lattice calculation. 

For the pseudoscalar, with the renormalization con- 
stant Zp = 0.21(2) determined in last section, the lattice 
calculation gives the result (0|P|0"+) « 8.6 ± 1.3 GeV 3 . 
It has been proposed that there is an approximate chi- 
ral symmetry between the scalar and pseudoscalar glue- 
balls sum rule is derived from an effective action 
which relates the topological susceptibility x m the pure 
gauge case to the gluon condensate Go, 

(60) 

where r\ ps 0.7. Using our lattice results, the degree 
of chiral symmetry rj can be obtained from the ra- 
tio of the pseudoscalar to scalar matrix elements [l6| . 
7] L = (0|P|0-+)/(2(0|S'|0 ++ )) re 0.3, which is also more 
than two times smaller than the result of QCD sum rule. 
These facts hint that there may be a substantial quench- 
ing effect in the matrix element of the scalar. 

We can also estimate the glueball contribution to the 
topological susceptibility by the lattice matrix element 
and glueball mass, 



Xglucball 



( (o|P|o- 



V327r 2 M -H 



(103 MeV) 4 ~ O.llx, 



( 61 ) 

which implies that the pseudoscalar glueball gives an ap- 
preciable 11% contribution to the topological susceptibil- 
ity. 
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TABLE XXI: The final glueball spectrum in physical units. 
In column 2, the first error is the statistical uncertainty com- 
ing from the continuum extrapolation, the second one is the 
1% uncertainty resulting from the approximate anisotropy. 
In column 3, the first error comes from the combined uncer- 
tainty of roMa, the second from the uncertainty of r$ = 
410(20) MeV 





r M G 


M G (MeV) 


0++ 


4.16(H)(4) 


1710(50) (80) 


2++ 


5.83(5)(6) 


2390(30)(120) 


0"+ 


6.25(6)(6) 


2560(35)(120) 


1+- 


7.27(4)(7) 


2980(30)(140) 


2-+ 


7.42(7) (7) 


3040(40)(150) 


3+" 


8.79(3)(9) 


3600(40)(170) 


3++ 


8.94(6) (9) 


3670(50)(180) 


1" 


9.34(4) (9) 


3830(40)(190) 


2"" 


9.77(4)(10) 


4010(45)(200) 


3" 


10.25(4)) (10) 


4200(45) (200) 


2+" 


10.32(7) (10) 


4230(50) (200) 


o+- 


11.66(7)(12) 


4780(60) (230) 



ball spectrum is summarized in Tab. IXXll and Fig. I1GI 

After the non-perturbative renormalization of the local 
gluonic operators, we finally get the matrix elements of 
scalar(s), pseudoscalar(p), and tensor operator (t) with 



12 



10 



+ 




2 + 




3 3 + 

















3 $ 

s 

CD 



In the tensor channel, the glueball matrix element is 
extrapolated to 1.0±0.2 GeV 3 in the continuum, which is 
the average of results of E and T2 channels. In the calcu- 
lation, it is found that in the lattice spacing range we use, 
the glueball mass and matrix elements are approximately 
independent of the lattice spacing, this implies that the 
lattice artifacts might be neglected here. If the renor- 
malization constant Zt ~ 0.52(15) of the tensor operator 
does not change much in the range of lattice spacing and 
applies to all the (3 values in this work, the renormalized 
matrix element of tensor operator is 0.52 ± 0.19 GeV 3 , 
which is in agreement with the prediction 0.35 GeV 3 from 
the tensor dominance model 0] and QCD sum rule 0] 
for the tensor mass around 2.2 GeV. 



VII. Conclusion 

The glueball mass spectrum and glueball-to-vacuum 
matrix elements are calculated on anisotropic lattices in 
this work. The calculations are carried out at five lattice 
spacings a s 's which range from 0.22 fm to 0.10 fm. Due 
to the implementation of the improved gauge action and 
improved gluonic local operators, the lattice artifacts are 
highly reduced. The finite volume effects are carefully 
studied with the result that they can be neglected on the 
lattices we used in this work. 

As to the glueball spectrum, we have carried out cal- 
culations similar to the previous work Q on much larger 
and finer lattices, so that the liability of the continuum 
limit extrapolation is reinforced. Our results of the glue- 



FIG. 16: The mass spectrum of glueballs in the pure SU(3) 
gauge theory. The masses are given both in terms of ro 
(r^ 1 = 410 MeV) and in GeV. The height of each colored 
box indicates the statistical uncertainty of the mass. 



the results 



s = 15.6 ±3.2 (GeV) 3 
p = 8.6 ± 1.3 (GeV) 3 
t = 0.52 ±0.19 (GeV) 3 



(62) 



where the errors of s and t come mainly from the errors 
of the renormalization constants Zg and Zt- The more 
precise calculation of Z$ and Zt will be carried out in 
later work. 
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